Nonergodisity of a time series obeying Levy statistics 
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Time-averaged autocorrelation functions of a dichotomous random process switching between 1 
and and governed by wide power law sojourn time distribution are studied. Such a process, called 
a Levy walk, describes dynamical behaviors of many physical systems, fluorescence intermittency 
of semiconductor nanocrystals under continuous laser illumination being one example. When the 
mean sojourn time diverges the process is non-ergodic. In that case, the time average autocorrelation 
function is not equal to the ensemble averaged autocorrelation function, instead it remains random 
even in the limit of long measurement time. Several approximations for the distribution of this ran- 
dom autocorrelation function are obtained for different parameter ranges, and favorably compared 
to Monte Carlo simulations. Nonergodicity of the power spectrum of the process is briefly discussed, 
and a nonstationary Wiener-Khintchine theorem, relating the correlation functions and the power 
spectrum is presented. The considered situation is in full contrast to the usual assumptions of 
ergodicity and stationarity. 
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I, INTRODUCTION 

Many time series exhibit a random behavior which can 
be represented by a two-state process 0]. In such pro- 
cesses the state of the system will jump between state on 
and state off. Examples include ion channel gating dy- 
namics in biological transport processes @, 0| and gene 
expression levels j3j_S in cells, neuronal spike trains Q, 
motion of bacteria [jj , fluorescence intermittency of single 
molecules 0] and nanocrystals [H El El El, and flu- 
orescence fluctuations of nanoparticles diffusing through 
a laser focus Some aspects of spin dynamics can 

also be characterized using two distinctive states 
These diverse systems may disp lay non-ergodicity and/or 
Levy statistics [ijj, [01 [01 l2i| , and often their behavior 
is found to deviate from simple scenarios used in the past 
to interpret the behavior of ensembles. In particular, in 
certain systems @, 0, 0, E3, E3, E3, E3, E3, EH E3| power 
law sojourn times are found for one or both of the states. 
Levy statistics, which manifests itself in appearance of 
power laws, is also found in flows on chaotic maps [2lj . 
which may be used to model dynamics of various com- 
plex systems with non-linear interactions. In this paper 
we address non-ergodicity of the Levy walk processes us- 
ing a stochastic approach. 

We model the intermittent behavior by a random pro- 
cess which switches between the two states after random 
sojourn times drawn from the probability density func- 
tions (PDFs) ip±(r), where the ± denote the two states 
(see Fig. QJ. It is assumed that these sojourn times are 
mutually independent random variables. In the follow- 
ing we assume common PDF for both states ip{r), unless 
stated otherwise, and assume that in state +, or on, the 
system is described by the intensity 1=1, while in state 
— , or off, it is described by zero intensity, I — (Fig. QJ. 
We consider the case of power law decay for long times 



t, 



time 



tp{r) - 6t~ 



< 6 < 1, 



(1) 



where we use natural units with dimensionless r. Such 



Figure 1: Schematic representation of a dichotomous process. 
T = T' — t' , where T' is the duration of the experiment and 
t' is the time difference used in correlation function (see Eq. 
J2J ). Note that in Section El we redefine t„ to be equal to T 
and r n is redefined to be T — t n -i, to simplify notation. 



distributions are observed in nanocrystal experiments 
[fll UnL ITlL ITiL IT^|. which under continuous laser illumi- 
nation exhibit random two-state blinking. As the mean 
sojourn time diverges, this situation reflects aging and 
non-ergodicity. Aging means dependence of some observ- 
ables (e.g., ensemble average correlation functions) on 
absolute times from the process onset at time zero, even 
in the limit of long times jH El El HE E3| Non- 
ergodicity means that ensemble averages are not equal to 
time averages of single realizations, even in the limit of 
long times. 

Generally speakingjOur model represents the so-called 
Levy walk process [19], in which a particle travels on 
a line with a constant velocity, changing directions at 
random times; the sojourn times are distributed with a 
power-law decaying PDF ip( T )- Some of the systems men- 
tioned above can in certain aspects be viewed as physical 
realizations of the Levy walk. 

In this manuscript we investigate the time average cor- 
relation function of the Levy walk process. When 9 < 1 
the process is nonergodic, because the mean sojourn time 
diverges. It is a common practice to replace the time 



2 



average correlation function with the ensemble average 
correlation function. Such a replacement is valid only for 
ergodic processes. Previous attempt to model correlation 
function of the Levy walk process, ignored the problem 
of ergodicity j2^|. Nonergodicity was observed in exper- 
iments of Dalian's group who obtained noner- 
godic correlation functions in experiments on nanocrys- 
tals. However, as far as we know there is no attempt to 
quantify the nonergodic properties of correlation func- 
tions of blinking nano-crystals and other Levy walk pro- 
cesses. Such a quantification is important in understand- 
ing the unusual behavior of physical systems and math- 
ematical models described in terms of Levy walks. Here 
we present a detailed analysis of our findings, part of 
which was reported in [29] . 



Q 
0- 



° <N)-100 

• (n> = iooo 
A iV r = 


: ° (N) = 100 
;< ■ <N) = 1000 

_4^Sy_ r - 0.01 


• <N) = 100 
(N>.1000 

JF r = - 1 


<N) = 100 
<N) = 1000 

, ~~s& jf , , ^ 


(N>.100 
(N> = 1000 


° (N) = 100 
<N) = 1000 



0.2 



0.4 0.6 

TA 



1 



0.2 



0.4 0.6 

TA 



II. TIME AVERAGE CORRELATION 
FUNCTIONS 

We consider an on- off signal in the interval (0, T') with 
intensity I(t) jumping between two states I(t) = 1 and 
I(t) = 0. At start of the measurement t = the process 
begins in state on 1(0) = 1. The process is character- 
ized based on the sequence {rf™, t%^, t™, , • • •} of 
on and off sojourn times or equivalently according to 
the dots on the time axis ii,i2, • • •, on which transitions 
from on to off or vice versa occur (cf. Fig. Define 
the following time-averaged (TA) correlation function for 
a single realization/trajectory: 



c T A(t',r 



i(t)i{t + t>)dt fii(t)i(t+t')dt 



and we denoted 



T> - V 



T = T - r > 0. 



T 



(2) 



We are interested in the asymptotic behavior of the cor- 
relation function for large T and t' , and define a ratio 

r = Yn (3) 

which will be a useful parameter. In the non-ergodic 
situations we consider, the distribution of the correlation 
function will asymptotically depend on t' and T" only 
through their ratio r. 

The mathematical goal of this paper is to investigate 
the PDF of C TA {t',T'). We first consider the PDF of 
CTA(t',T') in the ergodic case, and then address the 
non-ergodicity for 8 < 1. This PDF is denoted by 
Pc TA /t> ,t')( 2; )) where < z < 1 are possible values of 
C T A(t\T'), due to Eq. 0. 



A. Ergodic case 

Let us first consider the ergodic case with exponen- 
tial PDF of sojourn times ip(r) — e~ T , when the mean 



Figure 2: Distribution of time-averaged correlation function 
for ip( T ) = e_T is seen to approach the Dirac delta function as 
the average number of transitions per realization (N) grows. 
Location of the delta function shifts from 1/2 for r = to 
1/4 for any r / for large enough T'(and hence also t'), as 
indicated by the dotted line. Here (iV) = T' . 



sojourn time defined by 

/>oo 

(t) = / tiP(t)cIt = 1 
Jo 

is finite. If the process is ergodic, the PDF of Cta(^', T') 
will approach in the limit of long times T" — > oo, the 
Dirac delta function 



Pc TA{ t>,T>){z)~5{z-(C TA {t',T'))) 



(4) 



where (} represent ensemble average. This is what we 
mean by ergodicity of the two-time correlation function. 
We illustrate this behavior in Figure using numerical 
simulations. Increasing the experimental time T 1 (and 
hence also t', to keep r constant) leads to narrowing 
of the distribution of the correlation function, yielding 
asymptotically Eq. It is also clear that, for any 

nonzero r the ensemble average {CrA{t' ,T')) will tend 
to (1/2) 2 = 1/4 as we increase T'. Stretching of the dis- 
tributions observed in Fig. [2 for large r is due to the 
finiteness of T': here T = T' — t' becomes of the order 
of unity, which is the mean time of e~ T . Therefore, this 
behavior is completely pre-asymptotic. 

The picture is completely different when we consider 
Eq. (QJ with 9 < 1, as is shown below. There is no 
narrowing of the distribution, and it actually tends to 
a universal shape, which is a function of r and 9 alone. 
The analogue of this distribution in the ergodic case is 
the Dirac delta, Eq. (@J. In the ergodic case, one is 
usually interested in the non-universal behavior for rela- 
tively short t' of the order of mean sojourn time, while 
for long t' the behavior is trivial. On the contrary, in 
the non-ergodic regime we consider, the behavior of in- 
terest in this paper is the universal nontrivial asymptotic 
behavior. From now on, 9 < 1 |3Cj . 
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Figure 3: Two randomly selected trajectories for 9 = 0.3. 
There are approximately 1000 transitions in each trajectory. 
The behavior is dominated by a few large intervals and hence 
is strongly nonergodic. 




time 



Figure 4: One randomly selected trajectory for 9 — 0.8 with 
1000 transitions. In comparison to 9 — 0.3 (Fig. yj, the 
longest sojourn times here are shorter and the behavior is 
less nonergodic. 



B. Non-ergodic case 

We begin the discussion of a non-ergodic situation by 
illustrating two randomly selected trajectories for 9 = 0.3 
in Fig. UJ Clearly, these two trajectories are differ- 
ent, and hence time averaged correlation functions of 
these two trajectories will be different, yielding ergodic- 
ity breaking. It is important to emphasize that increas- 
ing the measurement time X", would not yield an ergodic 
behavior, since the process has no characteristic average 
time scale. In Fig. 0we show one trajectory with 9 = 0.8 
to compare to Fig. One can say that for 9 = 0.8 the 
nonergodicity is weaker. Unlike Fig. in Fig. 0]we do 
not see one long on or off period dominating the time 
series. In Figure El we plot ten typical realizations of a 
correlation function, for a power-law decaying ip(r) fol- 
lowing Eq. (JTJ with 9 = 0.3 and 9 = 0.8. The most 
striking feature of this figure is that the correlation func- 
tions are random. For very small r there is more or less 
smooth evolution of the correlation functions. As r grows 



Figure 5: Ten typical realizations of Cta dependence on r = 
t'/T' for 9 = 0.3 (top) and 9 = 0.8 (bottom). T' is kept 
constant, t changes. For an ergodic process all correlation 
functions would follow the same master curve, the ensemble 
average correlation function. 
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Figure 6: PDF of C T A(t',T') for different r = t'/T' and 9 = 
0.3. (N) w 10 3 , T' « 1.66 x 10 10 . Abscissas are possible 
values of Cta(£' ,T'). Diamonds are numerical simulations. 
Curves are analytical results without fitting: for r = Eq. 
is used (full line), for r = 0.01 and 0.1 Eq. (J2U is used 
(dashed) and for r = 0.5, 0.9 and 0.99 Eq. JSjJ is used (full). 
See Section IvTl for details. 



their behavior becomes more chaotic. We stress that this 
randomness is a true behavior and is not a problem in 
our simulations. 

For many realizations, our numerical simulations are 
used to obtain Pc T A (t< ,T')( Z ) depicted in Figures , 
and H for 9 = 0.3, = 0.5 and 9 = 0.8, respectively 
((N) is the average number of transitions per realization; 
details of these simulations are deferred until Section IVII 
and theoretical analysis is developed in SectionlVlbelow) . 
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The diamonds are numerical results. In all the figures 
we vary r = t'/T'. First consider the case r = 0. For 
9 = 0.3 and 9 — 0.5 we see from Figs. El and that 
the PDF Pc TA (t> ,T')( Z ) has a U shape. This is a strong 
non-ergodic behavior, since the PDF does not peak on 
the ensemble averaged value of the correlation function 
which is 1/2. On the other hand, when 9 = 0.8 the PDF 
Pc TA (t' ,T')( Z ) h as a W shape (cf. Fig. El, a weak non 
ergodic behavior. To understand the origin of this type of 
transition note that as 9 — » we expect the process to be 
in an on state or an off state for the whole duration of the 
measurement. This is so because the probability that the 
sojourn time is longer then X" will be ~ (T')~ e — > 1 (cf. 
Fig. EJ. Hence in that case the PDF of the correlation 
function will peak on C T a(*', X") = 1 and C T A(t' , T') = 
(i.e U shape behavior). On the other hand when 8 —* 1 
we expect a more ergodic behavior, since for 8 > 1 the 
mean on and off periods are finite, this manifests itself 
in a peak of the distribution function of Cta^ ',T') on 
the ensemble average value of 1/2 and a W shape PDF 
emerges (Fig r = 0). Note that for 8 < 1 there is 
still statistical weight for trajectories which are on or off 
for periods of the order of the measurement time T", and 
the distribution of Cta(0,T') attains its maximum on 
C TA (0,T') = 1 and C T a(0,V) = 0. 

For r > we observe in Figs. EJandElnon-symmetrical 
and non-trivial shapes of the PDF of the correlation func- 
tion. These PDFs agree very well with the analytical re- 
sults, which we derive later. Not shown in Figs. El Hand 
El is a delta function contribution on Cta{^ ,T') = 0. In 
other words, for t' ^ 0, some of the random correlation 
functions are equal zero. The number of such correlation 
functions is increasing when r is increased. When r — > 1, 
half of the correlation functions are equal to zero (see 
Section 0. Qualitatively, considering large r, the corre- 
lation is between the signal close to its starting point and 
the signal close to its end point. Roughly speaking, close 
to the end of the signal, typically long sojourn intervals 
with no transitions occur (cf. Fig. El i.e. persistence, as 
explained later in the paper in more detail - cf. Eq. (HJ). 
For those types of trajectories being in state off at the 
end, the correlation function should be zero. We stress 
that the distributions observed on Figs. El and El arc 
not a scaling artifact: analogous calculations in the case 
of 9 > 1 lead in the limit T" — > oo to Dirac <5-functions 
instead, as was shown above (Section |ll Al cf. Fig. EJ. 

We now turn to an analytical treatment of the de- 
scribed non-ergodicity. 



III. t' = 0: LAMPERTI DISTRIBUTION 
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Figure 7: PDF of C T A(t', T") for different r = t'/T' and = 
0.5. (TV) « 10 3 , T' « 2.47 x 10 6 . Diamonds are numerical 
simulations. Curves are analytical results without fitting: for 
r = Eq. is used (full line), for r = 0.01 and 0.1 Eq. O 
is used (dashed) and for r — 0.5, 0.9 and 0.99 Eq. 12911 is used 
(full). See Section E3 for details. 
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Figure 8: PDF of C TA (t',T') for different r = t'/T' and 6 = 
0.8. (N) w 10 4 , T' w 6.15 x 10 5 . Diamonds are numerical 
simulations. Curves are analytical results without fitting: for 
r = Eq. is used (full line), for r = 0.01, 0.1 and 0.5 Eq. 

is used (dashed) and for r = 0.9 and 0.99 Eq. J22j is 
used (full). See Section IvTl for details. 



the time average intensity between time a and time b > a. 
For t' = from Eq. J2J it immediately follows that the 
time averaged correlation function is identical to the time 
average intensity 



In the case t' = there exists known asymptotically 
exact expression for Pc ta (o,t) i z )- Let us define 



Cta{0, T) — T[o,t] = 



1 [0.T] 

T 



(6) 



•Ml 



b — a 



where Tj+ fe , is the total time on of a particular realization 

(5) 

in the time interval [a, b] . The time average intensity 
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1[o,T] has a known asymptotic distribution as T — > oo, 
found originally by Lamperti [113, lLU] and denoted in this 
paper as £$: 



Pc ta {o,t){z) =h (z)- 



and 



h (z) 



sin7r0 z - 1 {l-z) 



+ {1 - z) 2e + 2z (I - z) cos7t6 



(7) 

for < z < 1. For negative 2 and for z > 1 it is zero. 
Note that £e(z) = £e(l — z) and £g(z) diverges at z = 

0. 1. This function is normalized to 1 for any < 9 < 

1. The Lamperti PDF is shown in Figs. El □ and 
for the case r = 0, together with the numerical results. 
The transition between the U shape behavior and the 
W shape behavior happens at C = 0.5946.... Lamperti 
distribution is related to the well known arcsine law 
(case 9 = 1/2). Other works regardin g re lative time spent 
by a system in one of two states are flfiL l^lfl . \%<fy . 



IV. ENSEMBLE AVERAGE (C T A(t',T')) 

Another useful asymptotically exact result that can be 
derived is the mean of Pc TA (t' ,T')( Z ), i.e., the ensemble 
average of Cta(V ' ,T'). Generalizing to slightly different 
on and off time PDFs, with equal exponents but different 
coefficients, 



(8) 



it has been shown [26] that the mean intensity-intensity 
correlation function is asymptotically, for t' > 



(Ht)i(t + *)>~P+-P+P-—B 



1 



1 + t/t' 



1-9,9 



where the incomplete beta function is defined as 



(9) 



B(z;a,j3) = x a -\l - xf~ l dx (10) 
Jo 



and 



A± 



In the particular case of equal ip±(t) we have P± = 1/2. 
Eq. Q exhibits aging since the correlation function de- 
pends on t even when it is long. Aging of the ensemble 
average correlation function is related to nonergodicity 
of single realization trajectory. 

Integrating we thus obtain from Eq. ^ 



(C TA (t>,T>)) 



fi(I(t)I(t + t'))dt 



sin tt9 



T 



M- 



l-r 



(11) 



We see that the mean of the single trajectory correlation 
function asymptotically depends only on the ratio r of 
its arguments. We will show that the same is true also 
for the whole PDF of this random function, and not only 
for its mean. For r close to zero and to one, 



(C TA (t',T')) 



P + 1 - (1 - P + 



! sin tt9 



7T0(1 - 9) 



PI + Pi P_ 



(1 — r) e sin7r6 

7T0(1 +0) 



r < 1 



1 - r < 1. 



(12) 



It is worth mentioning that for an ergodic time series 
the variance 

a l( T ) = ((^[O.T] ~ (^[0,T])) 2 ) = (Z[0,T]) - { J [V,T\) 2 

should go to zero as T — > oo. In the case 9 < 1, in this 
limit (2[o,T]) - > P+ U$k an d using Eq. O, 

<4(T) -» S -^-P + P_ f 1 B(x; 0, l-9)dx = P+P_(l-0), 
t Jo 

(13) 

which is non-zero, and so we can prove the non- 
ergodicity of the considered process, even without know- 
ing Pc TA (t',T')(z)- The last equality can be easily ob- 
tained using Eq. 1101) . 

We conclude this section by introducing the probability 
Po(p>, b) of making no transition, either up to down or vice 
versa, between two arbitrary times a and b > a, known 
as the persistence probability. For large a (cf. Eq. l(B2j|) 

SI n TTr) 

p (a,b) B (a/b; 9,1-9). (14) 

Without going into details, we note that this probability 
plays important role in Levy walks, and in particular in 
formulas given above Its crucial feature is that it 

depends on the ratio of times and not on their difference, 
as is the case for ergodic processes. See also |3JJ,|3J|. 

Remark: Eq. l|T3"jl also follows from the fact that 
cr|(T) should approach the variance of the Lamperti dis- 
tribution (for P+ = P-), whose moments can be calcu- 
lated (lH appendix B]. 



V. t' / 0: APPROXIMATE SOLUTION 

We were able to obtain only a formal exact solution 
for the PDF of C TA (t',T') for t' ^ (see Appendix 
lA)l . Therefore, we resort to approximations. To start 
our analysis we divide the integration interval [0, T] into 
sojourn times Tj, For convenience we redefine the first 
tj > T to be equal to T, and denote its index by n: 
t n = T. Accordingly, t„ is redefined to be T — t n -i [cf. 
Fig. QJ], Thus, for i < n we write 



C T A{t',T') 



(15) 
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where we used the initial condition that I(t) = 1 at time 
t = 0. Hence I(t) — 1 in <t <U when i is odd, oth- 
erwise it is zero. The summation in Eq. l(T5|l is over odd 
i's, and t„ = T, namely n — 1 in Eq. 11511 is the random 
number of transitions in the interval [0, T], From Eq. 
itlfil) we see that the time averaged correlation function, 
multiplied by T, is a sum of the random variables 

U ( Ti - t' +T^. tt . +t >]t' i odd, n > t' 

I(t)I{t+t')dt = I T {t ._ l+VM+ t')n i odd, n < i! 

I i even. 

(16) 

Using Eqs. itlfil ITfil) we find an exact expression for the 
correlation function 

n n 

TC TA (t\T')= J2 (i-W,^)^ 

i odd j odd 

r* < f 

n 

- *' J2 i 1 -Ziu.u+f])- 

i odd 
n > t' 

(17) 

The first term on the right hand side of this equation is 
T + the total time spent in state on in the time interval 
[0,T], in the remaining two terms we have considered 
sojourn times Tj larger or smaller than t' separately. 

The core idea of our approximate solution is to replace 
the time-averaged intensities entering Eq. itTTIl by their 
mean-field value, specific for a given realization. Then for 
short H we replace I [U l+t , u+tl] and I[ ti , ti +t'] b Y Z[o,T], 
while for long t' we use 1\t> t T'] instead. Some alternative 
approximations are given in Appendix [O In the follow- 
ing, we treat short and long t' separately. 



A. Small t' 



Within the mean field theory, Eq. ifTTjl is approximated 

by 

TC TA {t\ T') = 1 %T ]T- (1 - X [0 ,t]) {t'N+ + E+) (18) 

where N + is the number of odd (i.e. on) intervals satis- 
fying n > t' and i < n, while S + = Z)" dd,r ( <t' Ti is the 
sum of all odd Ti < t' and i < n. For any particular re- 
alization N + will decrease with t' in a step-wise fashion, 
while S + will increase in a step-wise fashion. The term 
t'N + + S + in Eq. I|18J1 . however, will be continuous. 

We proceed by replacing N + and S + with their scaling 
forms. N + should scale as ~ n + J t , tP(t)cIt and E + ~ 

n + J Q tiP{t)oIt, where n + is the number of on intervals 
comprising a given T + . First note that for t' > T + , 
N + = and E + = T + . Second, we assume 



sin 7Tfc 
tt6» 



(T+) 6 



(19) 



in analogy to the scaling of n with T (e.g., [la])- There- 
fore, using Eq. Q we propose that for 1 <C t' < T + 



N+ 



sm Tit 



and similarly, 



- 1 



(20) 



(21) 



Finally, plugging Eqs. into Eq. itT^ll results in 



C T A{t',T') 



Mo.t] 
I 2 

[0,T] 



L [0,T], 



t 1 r )^-[o,T] 



it8 (l-r)I [0T] 



t' < ^^ 
t' > ^^ 



(22) 



Eq. <22l) yields the correlation function, however unlike 
standard ergodic theories the correlation function here is 
a random function since it depends on 1[o,t] ■ 

The PDF of C T A(f ', T') = z is now easy to find from 
the Lamperti PDF of 2T[o,t] = x. Using the chain rule, 
and Eqs. <EE|[22j): 



Pc T A(t',T')(z(x)) 



dz(x) 



dx 



(23) 



which is a parametric representation of Pc TA (t' ,T'){ Z ) 
(dz/dx = dC T A(t',T')/dl l0 ,T] is found from Eq. ffi)t ). 



In Figs. H □ and we plot the PDF of C T A{t\T') 
(dashed curves) together with numerical simulations (di- 
amonds) and find excellent agreement between theory 
and simulation, for the cases where our approximations 
are expected to hold r < 1/2. In the above treatment we 
approximated I[ t ' ,T'} by 2T[o,t]i which is legitimate only 
for small enough t' < T, leading to a deterministic de- 
pendence of CrAit' ,T') on 2T[o,t]- 

Remark 1: Note that in the ergodic case (in which 
we can insert 6 = 1 in the scaling relations) it fol- 



lows that Cta(0,T') = J, 
Cta(0,T') 



T 2 

[0,T] 



[o,T] = 1/2 for r = and 
1/4 for any r ^ 0. This behav- 
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ior reflects complete decorrelation of I(t) and I(t + t') for 
any (large enough) t', irrespective of the value of T', as 
is indeed the case. 

Remark 2: There is a certain similarity between Eq. 
Il22t and Eq. 112t for small r. Only qualitatively, a re- 
alization with a given 2j ,T] can be viewed as generated 
using iP±(t) with A + A- (cf. Eq. JHJ), such that 
P+ = I\q,t\ ■ See additional discussion of Eq. H22II in 
Appendix IdI 



B. Large t' 

To understand the behavior of the PDF of the cor- 
relation function for the limiting case t' « T" 3> T the 
concept of persistence is important (see Eq. 114|l ). Recall 
that the probability of I(t + t') — const on the interval 
[t',T'] grows to unity as t'/T' — > 1. Moreover, there is 
virtually no dependence on the signal values on t G [0, T] 
and thus 

ifcxACf ,T')(*) « + - 0). (24) 

There is a collapse of half of the trajectories to a <5-peak 
at z = 0, because of zero intensity of the signal on [t', T'] 



in one of the two states, with probability — > 1/2. In 
the second case the signal will be unity throughout the 
interval [t',T'], with probability — > 1/2, while its relative 
on time distribution in [0,T] is given by Lamperti PDF. 

More generally, for t' not so large, but still t' > T we 
use the mean-field, or decoupling approximation yielding 
from Eq. itTTjl 



C T A(t',T') 



' Z[0,T]Z[t',T']- 



(25) 



To calculate the PDF of C TA {t' ,T') in Eq. fiB we 
use two steps: (i) calculate the PDF of Tw t 1 ] = z 
which statistically depends on Z[q,t] (it is denoted as 
Piu, T /, { z \1[q,t})) and then (ii) using the distribution of 
^[o,T], which is the Lamperti's PDF Eq. JJJ, calculate 
the PDF of Cta(^', T") = z: 



Cta( 



(v)W~ / ^(^[^,(^1*) 



fi;j; 



(26) 



Using the persistence probability Eq. 111411 . we ap- 
proximate the conditional PDF of I\t',T'} = z f° r a given 
1[o,t] m the case T <C t' by 



P V , T ,, (z|X[o,T]) ^ [1 - po (T, T')] Qx [t ,, T ,, (z) + p (T, T') [X [0 , T] 5 (z - 1) + (1 - J %T] ) 6 (z)] , (27) 

where Qi [t , T ,, (z) is the PDF of Z[t',T'] conditioned that at least one transition occurs in [T,T'], In Eq. I|27jl we 
introduced the correlation between I[t',T'] and Z[o,t] through the dependence of the right hand side of the equation 
on I[o t T]- We assumed that in the case of no transitions in the time interval [T,T'], the probability of the interval 
\t! , T'] to be all the time either on or off (the only possible choices) is linearly proportional to the value of I[o,t]- 
The persistence probability controls also the behavior of 

Qx M (z) ~ [1 -po (<', T')} 6 (0 < z < 1) + Pa (f, T') 6 (z) + S - (z ~ 1} . (28) 

We assumed that if a transition occurs in the interval [f, T'] the distribution ofTw t T'] is uniform [i.e., 6 (0 < z < 1) = 1 
if the condition in the parenthesis is correct]. This is a crude approximation which is, however, reasonable for our 
purposes (however when 9 approaches 1, this approximation does not work). The delta functions in Eq. II28|) arise 
from two types of trajectories: If no transition occurs either I[{',x'] — 1 (state on) o r Iw ,T'] — (state off) with 
equal probability. An asymptotically exact expression for Qi , T , (z) is given by Eq. l|B3t in Appendix El given the 
approximate nature of our derivations, however, we chose to use Eq. 12811 because it is much simpler. 
Finally, from Eqs. H27I28I26|) . and using S(a/x) = x5(a) for x > 0, we obtain after some algebra 



b M ( V ) (*) — [1 — Po (T, T')] ( [1 - p (t>, T')] J] l -^dx + P -^± [l g (z) + S (z)} \+p (T, T') \zl 6 (z) + 



(291 



Note that to derive Eq. 112911 we used the fact that 1[o,t\ these approximations are expected to hold, r > 1/2. Eq. 
and 2jt',T'l are correlated. In Figs. 16171 and |U we plot 12411 is recovered from Eq. l|29Jl in the limit of t'/T' — » 1. 
these PDFs of CrA(t',T') (solid curves) together with 
numerical simulations (diamonds) and find good agree- 
ment between theory and simulation, for the cases where 
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VI. NUMERICAL SIMULATIONS AND 
COMPARISON TO APPROXIMATIONS 

We performed Monte Carlo simulations to generate 
distributions of the time averaged correlation function 
Cta^' ■, T') for different values of r — t' /T' and with dif- 
ferent 9. Specifically, for each chosen 9 the function 



T > 1 



0, T < 1 



was used to generate random sojourn times until certain 
cumulative time T" 3> 1. This constitutes a single real- 
ization. Tens of thousands of realizations were generated 
for each 9. 

For each realization, Cta^ ', T') was calculated for dif- 
ferent t' using Eqs. <j2J and l|A2|l . To check whether the 
PDF of CTA(t',T') depends only on r we used different 
T' . We also used the one-sided Levy PDF for if){r) and 
found that our results do not depend on details of ip(r) 
besides the exponent 9 of course. In addition, we calcu- 
lated (CTA(f ,T')) from our simulations and compared 
it to the theoretical result Eq. l(TT|l . The agreement is 
excellent, as long as T = T' — t' 3> 1. 

Some simulations are shown on Figures and |H| to- 
gether with various theoretical approximations, for 9 = 
0.3,0.5 and 0.8 respectively. (N) is the number of tran- 
sitions made until time T", averaged over realizations. 
Diamonds are simulated data. Solid lines for r — are 
£e(z) where < z < 1 are possible values of Cta[0, T'). 
Dashed and solid lines for r / are Eqs. l(23jl and l(29j) 
for r < 0.5 and r > 0.5, respectively. 

The discontinuity of the dashed lines, which can be 
noticed at small values of Cta^ \ T') for r = 0.1 is due to 
the discontinuity of the derivative in Eq. (1231 at I[o,ti = 
r/(l — r), when Cta{V \T') becomes equal to T? Q r , = 



r 2 /(l 



', which is very small for small r. Overall, 



however, Eq. lf23j) agrees with the shown simulations for 
r < 0.5. 

Approximation ll29l) works well for all 9 values and 
r > 0.5, for which it was designed, and it can be seen 
that as r grows toward 1, the asymptotic result Eq. Jl2 II) 
is approached. The assumption of uniform distribution 
of 1[t>,T'] f° r values between and 1, used in Eq. p8f . 
is an oversimplification when 9 = 0.8, which is partly re- 
sponsible for slight discrepancies with the simulated data. 
Qualitatively, the PDF of I[ tl t2 ] is similar to £g(z) which 
starts growing a maximum at z — 0.5 for approximately 
9 > 0.6 and so Eq. l|28|l is not very accurate for 9 = 0.8. 
Also, here T'k6x 10 5 is not very large and therefore the 
simulated distributions haven't completely reached their 
asymptotic forms (e.g., observe slight shape differences 
between simulated data and theory for r = 0). 

Dot-dashed lines in Fig. E|for 9 = 0.3 are based on Eq. 
l|C3jl . They are shown only for r < 0.5; this approxima- 
tion works well in the limit of small 9 and r. 

Our simulations show that for 9 = 0.5 the 
PDF of Cta(1' ,T') = z is closely approximated by 



2 (C TA p, T')) 4.5 W + - 2 (C TA (t', T')))5{z). Dotted 
lines in Fig. are the nonsingular part of this expres- 
sion, i.e., Lamperti distributions normalized by the rela- 
tive mean. They are in good agreement with the data, 
and therefore are hardly visible. We have no explanation 
for this fact, besides the qualitative argument that as r 
grows from zero, for lower 9 the left side of the distribu- 
tion drops (cf. Fig. EJ, while for higher 9 it rises (cf. Fig. 
IHJ, and so somewhere between 9 = 0.3 and 9 = 0.8 it 
might remain unchanged. For t'/T' — » 1 this expression 
approaches Eq. J23J. 

For r = the PDF of C T A(t',T') is the Lamperti dis- 
tribution. As can be observed from comparison of the 
PDFs with r = and r > in Figs. EJandU the PDF of 
CTA(t',T') is shifted to the left as r increases from zero. 
This is so because small Z[o,t] values mean small propor- 
tion of time spent on, and there is a large probability 
that a small i! will yield zero correlation in such realiza- 
tions. This is in agreement with Eq. l|C3|l . For larger 9 
values, there are more short and less long intervals cov- 
ering the time of "experiment" T" (as can be seen from 
Eq. l|C2|t . because r 1_e increases toward 1 with growing 
9, for any fixed r > 0). Therefore, relatively small shift 
t' (small r) will cause no significant effect in the case of 
small 9, dominated by large intervals, while in the case 
of large 9 this small shift t' will decorrelate many inter- 
vals, thus significantly reducing the correlation function. 
Realizations with small 2"[q,t] also will lose correlation 
faster for the same reason, leading to a non-uniform vis- 
ible deformation of the shape of the PDF of Cta($ ', T'). 
Of course, as r grows this simple picture breaks. How- 
ever, for r approaching unity we recover another simple 
asymptotic result (124L 



A. 2D histograms 

Two dimensional histograms, showing the frequency 
of events CTA(t',T') for a particular value of I[o,t'] are 
now considered. These histograms show the correlation 
between Cta^' ■, T') and I^t']- As we explained already 
for r = we have CTA{t',T r ) = I[o,t']> hence we have 
total correlation in this simple case. When r is small, 
our approximate solution Eq. Il22t suggests a strong cor- 
relation between Cta^' ,T') and 1\q,t '] < However, the 
arguments we used to derive Eq. j23 neglect fluctua- 
tions since they are based on our non-ergodic mean field 
approximation. To check our mean field, and to under- 
stand its limitations, the two dimensional histograms we 
consider in this section are very useful. In addition, for 
large r we see from Eq. l(29jl , that according to the decou- 
pling approximation, the correlation between Cta(^', T') 
and T[o,t'] ls expected to be weak, as is demonstrated 
indeed by correlation plots in Fig. El 

Leaving the details to Appendix |E1 we can derive the 
following rigorous boundaries (i.e., the inf and the sup) 
of C T A{f,T'): for r > 1/2 
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Figure 9: Distribution of Cta(^! ,T') as a function of 2"[o,t'] for different values of r and 6. The gray scale is changed 
logarithmically with the number of occurrences inside a square bin. Darker regions mean higher occurrences. Dashed lines are 
Eq. 122H with I[o,t'] used instead of 2"[o,t]- Full lines CrA(t' ,T') — (2T[o,t'] ~ r) /(l — r) are shown as well. 



max < 0, 



X, 



[0,T'] 



1 



<C T A(t',T') <minh, 



L [0,T'] 



2(1 



r) 

(30) 

which is in agreement with Fig. El For 1/3 < r < 1/2 we 
obtain 



max < 0, 



L [Q,T'] 



hr-1 21, 



1-r 
< min 



<C TA (t\T') 



■ f 2 1-10,T>] ^[O.T']+ 1 — 2r 1 

m \3(l-r)' 2(l-r) / 



(31) 



For small r, on one hand the matters become more 
complicated, so our argument is more qualitative . First 
notice that if TZ T „ = T' — Xj+ T ,j is small enough then 
all the off intervals can lie inside [t 1 ', T] and be used twice 
(once in I(t) and once in I(t + 1'), for r ^ 0) to multiply 



r). Compare this to Eq. l|C3Jl derived in Appendix [O 
It is argued there that this value of CTA{t',T r ) will be 
achieved more often for lower 9, in agreement with Fig. 
El This is not a rigorous upper bound, though; see Ap- 
pendix|El If T^ T ,^ is small enough then the lower bound 
will be zero. The sufficient (but not necessary, in general) 
condition to achieve zero is T^ T n < 1/2 as then we can 
construct a trajectory by choosing zero intensity at the 
time t + t' if it is 1 at time t, and vice versa. 

On the other hand, for very small r (but t' can be 
large) we know that Cta^^T 1 ) is almost unchanged, as 
the whole signal is dominated by relatively few largest 
sojourn intervals (cf. Eq. gSJ); hence C TA (t',T') will 
be close to 2[o,t']- The rigorous bounds are, therefore, 
hardly reached. 



Remark: Our approximation Eq. i!22fl . with X\ 



L [0,T] 



on intervals, hence C TA {t' ,T') > (t - 2T [Q r ,J /T = replaced by T[o,T'\, is shown by the dashed lines on Fig. 



[2T[ T i] — r — l) /(l — r). In most cases, small T, 



[0,T'] 



means that the last sojourn interval (going up to time T") 
is in state on and all the off intervals are inside [0, T], so 

that C TA (t', T<) <(T- Tj-^) IT = (I [0 ,t>] - r) /{I - 



El For r = 1/2 it actually reduces to C TA {t', T 1 ) = lf TI] , 
which works better for higher 9, when the non-ergodicity 
is weaker. In fact, a more precise way to find Cta^' , T') 
in this case is using Eq. il25|l . but then there is no simple 
formula connecting Cta^' ',T') and 2"[q,t] like Eq. II22H . 
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B. Power spectrum 

It is useful to look at power spectra (PS) of generated 
intensity signals [IE H3, H3, HH ■ Power spectrum is de- 
fined as 



S(lj) = 



/M/(-") 

T' 



(32) 



where I(u>) is Fourier transform of I(t) (cf. Eq. f Afif 1 . 
We calculate such PS and find, as expected, that they too 
exhibit a nonergodic behavior, as shown in Fig. Each 
PS is random and does not fall on the ensemble averaged 
curve (dashed line) even after averaging the data in large 
frequency windows. Note that for smaller 9 the PS values 
for a given uj are spread wider, which is a reflection of a 
wider distributions of correlation functions (cf. Figs. El 
EJ). In light of the scaling C T A{t',T') ~ A - Br 1 - in 
expressions ltl2l and l(22)l for small enough r (but for t' 
as large as desired, as long as T' is large enough) we can 
argue that the PS will scale as (cf. Eq. &AF 



S(u) 2B(T') e - 1 Re/ T (t') 1 ~ e e~ iuJt ' dt' 

w 2BT' cos(7r6»/2)r(2 - 9)(luT') 9 - 2 oc uj 9 - 2 

(33) 

as long as w > 1/T" (term A in CTA(t',T') leads to a 
term 2 AT' sin wT'/(wT') which is zero for all u ^ used 
in calculating discrete power spectrum). This is indeed 
the case, as illustrated in Fig. EH In Eq. (l33f) we esti- 
mated the PS by Fourier transforming the correlation 
function, implying the well-known Wiener-Khintchine 
theorem. This theorem, however, is assumed valid only 
for stationary processes and for ensemble averaged cor- 
relation functions and spectra. Nevertheless, one can say 
that with respect to short sojourn times each realization 
is identical, and the observed non-ergodicity is due to 
necessarily poor statistics of long intervals (leading, in 
particular, to different values of B for different realiza- 
tions). See Appendix El for discussion on a generalized 
Wiener-Khintchine theorem. 

For the ensemble-averaged spectrum, using the value 
B = sin(7T0)/(47T0(l - 6>)) from Eq. Ip). 



(S(uj)) cos(ir9/2) 



(ujT 1 ) 9 - 2 , wT'»1. 



(34) 



T 2r(l + 9) 
This line is also shown in Fig. 

VII. SUMMARY 

We investigated autocorrelation of a dichotomous ran- 
dom process governed by identical waiting time distribu- 
tions of its two states, characterized by zero and nonzero 
intensity. We considered the case of a power law wait- 
ing time with exponent 9 lying between and 1, as this 
choice is of considerable practical interest. This process 
is a one-dimensional Levy walk process. Such power law 
distributions are experimentally observed, as discussed 




10 10 
1+coT 



Figure 10: Power spectrum S(u>)/T' for ten typical realiza- 
tions shown in Fig. |&] for 9 — 0.3 (top) and 8 — 0.8 (bottom). 
Data for each realization are averaged in exponentially in- 
creasing with uj bins. Each curve is normalized in such a way 
that at uj — the PS equals lf 0T ,y The PS is random due to 
non-ergodicity of underlying process. Dashed lines are given 
by Eq. 134ft and scale as ui e ~ 2 . The abscissas are 1 + luT' in 
order to show the value of PS at zero frequency on a log-log 
plot. 



in the Introduction. These distributions lead to aging 
and non-ergodicity and in particular, to a distribution 
of possible values of a single trajectory two-time corre- 
lation function for fixed times, even in the limit when 
these times go to infinity. This is in striking contrast 
to the standard situation in which correlation function 
asymptotically assumes only one possible value for fixed 
times, equal to the ensemble average (ergodicity). 

For our theoretical analysis of distributions of correla- 
tion functions we used the non-ergodic mean-field and the 
decoupling approximation, Eqs. iflSl) and m which 
various temporal averages of the intensity were replaced 
by the total time averaged intensity I[o.t] or 2[t',T'l) spe- 
cific for each realization. We then expressed the correla- 
tion function as a (deterministic or random) function of 
this time average. This enabled us to derive approximate 
results for the distributions of correlation functions from 
known distributions of time averaged intensity. We also 
related power spectra of single trajectories to the time 
averaged correlation functions, and demonstrated their 
nonergodicity as well as universal scaling which is a func- 
tion of the exponent 9 only. Our results agree well with 
numerical simulations, and, importantly, clarify the na- 
ture of the investigated non-ergodicity. Generalizations 
of our approach to situations with different on and off 
time distributions are possible. 

In the context of blinking nanocrystals, we showed 
HE] that the exponent 6 — 1/2 is a result of a simple 
model of first passage time of charge carrier in three di- 
mensions, based on standard diffusion. The experiments 
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[lfl H^. FLU lift l4ll | show, that rather generally, power 
law sojourn times describe dynamics of single particles 
in diverse systems. Since power law sojourn times (not 
necessarily for a two state process) lead to non-ergodic 
behavior, we expect that stochastic theories of ergodicity 
breaking will play an increasingly important role in the 
analysis of single particle experiments. 
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Appendix A: FORMAL SOLUTION 

We express the numerator in Eq. (2) through the 
cumulative renewal (transition) times tj by noting that 
I(t) = const on intervals t n < t < t n+ i, while I(t + 1') = 
const on intervals t m —t'<t< t m +i — t' and keeping in 
mind the restriction < t < T. For convenience, we rede- 
fine the first tj which is > T' to be equal to T". The index 
of this tj is denoted by N. Note that tn = T' — ijv-i is 
not distributed according to The temporal dura- 

tions (lengths) of intervals where both I(t) and I(t + t') 
are constant, are then 

l mn (t') = max {0, min(i n+ i, t m +i - t') - max(t„, t m - t')} 

(Al) 

and in particular 

l nn = max{0,T„ + i -t'}. 
Obviously, l m <n,n = and hence 

„T'-t' N-l N-l 

/ I(t)I(t + t')dt= Y] V l mn . (A2) 

n = I), m = n, 
n even m even 

Here, and throughout the article, we assume that the 
process starts in state on. This assumption is clearly 
asymptotically negligible, and is made here simply for 
purposes of notation. 

Using the cumulative PDF of {ri, tn-i I and N — 1 
under the constraint ijv-i < T' 15j (and because 
Prob[<Ar_i = T'] = 0) we can write formally the PDF 
oiTC T A{t',T') = x as 



/JV-l 



TC TA (t',T>) 



(,:) 



Hrk)dr k V(t N = T')x 



N-l 



N-l 



xG(T'-tjv-i)«5 



^ t y t lmn(t ) 

n = 0, m = n, 

n even m even J 



where 8 is Heaviside step function, V{tN — T') is the 
probability that no transition occurred between t^-i and 
T', and 6 is the Dirac delta. 

In order to get rid of the max and min functions in Eq. 
(I A 111 , one can perform Laplace transform of Eq. (I A If) 
with respect to t! (f — » u) and write similar expression 
for the PDF of TCta(u, T') (for real u). It can be shown 
that 



^m>n,n (^) 



and 



Inn (^) 



Tn+1 



-UT n+1 



(A3) 



(A4) 



We could not, unfortunately, utilize these expressions to 
calculate fTC TA (t',T')[x) or f T c TA (u,T>)( x ) and therefore 
have to resort to various approximations. 



Relation to power spectrum 

We derive a generalized form of Wiener-Khinchine the- 
orem for nonergodic nonstationary processes. In analogy 
to the numerical spectral analysis of a time series, we 
assume here that the intensity signal is identically zero 
outside of the interval [0, T']. Then the Fourier transform 
of intensity is defined as 

I{t)e- luJt dt= I I{t)e- Wt dt (A5) 

-oo J0 

and the power spectrum of a realization is defined in Eq. 
P2f. From Eqs. (SB and © 

pT' pT' 

T'SV) = / I(ti)e- iutl dt! / I{t 2 )e iut2 dt2. 
Jo Jo 

We now divide the integration over t 2 into two parts and 
replace the order of integration in the first part: 



T' pT' pT' pt ± pT' pT' 

dti I dt 2 = / dti I dt 2 + dti / dt 2 
i) Jo Jo Jo Jo Ju 



T' pT' pT' pT' 

dt 2 / dt x + dti \ dt 2 . 
io Jt 2 Jo Jti 

Swapping names t\ and t 2 in the first part thus yields 

T'S{lu)= ( dti [ dt 2 I(t 1 )I(t 2 )[e lu( - tl - t2) +e t ^ t2 ~ tl ' > } 
Jo Jti 

and introducing t = t\ and t' = t 2 — t\ results in 

r T' rT'-t 



1 



o ./o 



l\[~— iiijt' i Aiot 1 ~\ 



S{lu) = — I dt I dt'I(t)I(t + t') [e~ wt + e Mt ] 
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In a similar fashion, we can write the Laplace transform 
of 

K(t , ,T') = (T'-t , )C T A(t',T')= f ' dtl(t)l(t + t') 

Jo 

(A6) 

with respect to t' as 
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K(u,T') = I dt'e~ ut I dtl(t)l(t + t') 
Jo Jo 



T'-t' 



t' /.r'-t 

dt \ dt'I{t)I(t + t')e 
Jo 



l\ n -ut' 



and it becomes evident that 

T'S{u) =K(iw,T') + K(-ioj,T'). (A7) 

This is a generalization of the Wiener-Khintchine theo- 
rem stating that the power spectrum is given by cosine 
Fourier transform of a correlation function. But while 
this theorem is used for ensemble-averaged correlation 
functions of stationary processes, here we have a sim- 
ilar relation for a non-stationary process, and without 
ensemble averaging. Note that the dependence on T' is 
preserved, in contrast to the regular Wiener-Khintchine 
theorem. 

From Eq. ijAfijl it follows that 

v( r p , \ vis, ( rp/N . dC T A(u,T') 

K(u, T ) = T Cta{u, T ) H 

ou 

and for large uT > 1, C TA (u,T') « K(u,T')/T', lead- 
ing finally to 

5(w) = C T A(iu>,T') + CTA(-iu,T') 



d[C TA (iu,T') - C TA (-iu,r)] 
idu 

C T A(iuJ,T') + C T A{-iu,T') 



(A8) 



for large vT' 3> 1. Note that for a single trajectory 

correlation defined as ( J Q T * I(t)I(t + t')dtj /T' instead 

of Eq. (|2j| the generalized Wiener-Khintchine relation is 
exact for any uj (cf. Eq. O). 

As an illustration, consider now our case of the on- 
off process. Fourier transform of an intensity I{t) for a 
realization is: 

N-l 

I(u) = V e - iutn+1 (l - e^+M. 

n = 0, 
n even 

Then it is straightforward to show that 

T'S(lu) = K(iu, T') + K{-iuj, T'), 

wher e K (u,T') can be found utilizing Eqs. \A2l . (|A3f) 
and l!A4|l . This is a particular case of the general relation 
ifATIl . 



a, 6] 



Here we present asymptotically exact formula for a dis- 
tribution of on times on an arbitrary interval [a, b], where 
a and b — a are large enough. We denote the first renewal 
time after a by v. We have to take two possibilities into 
account. First is that there was at least one renewal in- 
side the interval and then a < v < b. Thus 



1 [a,b] 



Y 



where Y is the on time from a till first renewal v. Asymp- 
totically, Y is independent of initial conditions and its 
PDF is 

fY{y) = \6{y-{v-a))+ 1 -8{y-Q), 

where the two Dirac deltas correspond to being in state 
on or off. 

After renewal at time v, again asymptotically, we can 
use the PDF of which is also independent of its 

"initial condition", i.e., the value of Y being or 1. Then 
the PDF of Tj+ b J(b— v) is given by the Lamperti ie and 
therefore for any fixed v 



f T + (x\v;v<b) 



2(6-i/) 



(Bl) 



The second possibility is that v > b and in this case, 
clearly, 



f T + (x\v\ v > b) 



l -5(x-{b-a))+ l -8{x-0). 



Introducing the PDF of the forward recurrence time, 
fs{v — a; a), which is the PDF of having to wait for the 
first renewal after time a for a period of time v — a [THI ] . 
we finally obtain 

f T+ (x) = f T + (x\v)f E {v - a;a)dv 

K&] J a la,b] 

POO 

+i (S(x - (b - a)) + 5(x - 0)) / f E {v - a; a)du, 



with f T + (x\v) in the first integral given by Eq. IB II . 

[a, fa] 

The last integral 



Po(a,b) 



fftiy — a; a)dv 



(B2) 



defines the persistence probability po (a, b), for b> a. The 
function $e(v — a; a) is equal to i>(v) for a — 0. For large 
a, it is the Dynkin function (lH Eftl32l | 



f E (f - a; a) 



sin tt9 



and then po(a, b) can be written as in Eq. 11411 . 
Finally, the PDF of Z[ a ,&] = z is 

Q I[ab] {z) = {b-a)f T+ ((b-a)z). (B3) 

[a,b] 
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Appendix C: PARTICULAR SOLUTIONS 

We consider here two situations which can be analyzed 
differently from the approach presented in Section 
This analysis helps understanding the structure of the 
correlation functions. 



1. Extremely small t' 

If t' < t%, i = 1,...,N then, using Appendix El we 
obtain l m>n +i, n = 0, l nn = T n+i - t' and 



the "experimental" time T . Then, for m > n using 

m 



m+l 



tn+1 — (tm+1 — t') — f — 



i=n+2 



TC TA (t',T') = T, 



[o,t>; 



N+l 



, Nt' 
1 ~ 1 [a.T'] 2 

(CI) 



yields 



also 



1+0 



Inn — dn,k-l(T~k — t'), 



where the subscript int indicates that [...] denote integer where is Kronecker delta, and hence 
part. It is easy to see that in this case, T£ T ,, > for 



" [0,T'] 

t' > and, moreover, < TC T A{t',T r ) < T, as it should 
(because N > 1). 

The fraction of time T" covered by short intervals < 
t' scales as 



Jo' Mt)dt 



(C2) 



for large enough i' and T'. Hence the contribution of 
these short intervals is negligible if t' <C T", although i' 
is large (in contrast to the case when the mean sojourn 
time is finite and the fraction of time covered by intervals 
shorter than t' grows to 1 as t' increases, irrespective of 
the ratio t'/T'). Therefore, we argue that Eq. IIClJl can 
still be used when t' T', if by N we understand the 
number N e ff of intervals longer than t' . It is important, 
however, to distinguish these coarsened intervals from the 
original intervals T{ . The durations r e f f of the coarsened 
intervals are not governed by ip(r). For a power law V , ( r ) 
we expect that their durations are still governed by a 
power law PDF with the same exponent 9. Nevertheless, 
it is questionable to use asymptotic expressions for N e f f 
as a function of t'/T', constructed in a fashion used in 
Section because the PDF of r e // does not have to be 
a power law for r e // ~ t' (and it is zero for r e // < t'). 



2. Small 6 and intermediate t' 

It is possible to make an exact calculation if there exists 
an interval number k, and for t' such that 

N 

n<t' < T fc . 

i = l 
i ^ k 

Ubiquitous realization of this condition could be expected 
for small 0, when the longest interval often approaches 



Cta^X 



0. 



L [0,T' 



1-r 



k even 



-, k odd. 



(C3) 



In the case of odd k, Tj+ T ,j > r^. > t' so that always < 

Cta(P' , T') < 1, as it should. This particular solution 
also plays an important role in defining the boundaries of 
the two-dimensional correlation plots discussed in Section 

ED 



Appendix D: NOTES ABOUT EQ. (f22l 

In this appendix, we discuss some approximations in- 
volved in the derivation of Eq. l|22ll and some of its short- 
comings. 

We begin with Eq. 119(1 . Scaling behavior n oc T 9 , 
where n is the number of transitions up to time T, is 
well-known for < 9 < 1 (e.g., However, the 

distribution of n is wide and its standard deviation is 
also known to scale as T e . For our purposes, we want 
to represent this standard deviation as arising from two 
contributions. First is that n depends on T + = Tj+ T j, 

while second contribution is that for any fixed T + there 
still is a distribution of n values. We can approximate 
the first contribution by writing n oc (T+(T -T + )/T)°. 
Since T + cx T, this formula does not contradict standard 
scaling n ocT e , and it is at least in qualitative agreement 
with our numerical simulations. To justify it we observe 
that when T+ T then there is probably a large in- 
terval of state off, which covers almost all the time T. 
If we remove this large interval then the remaining total 
time will be of the order of T + , while the number of in- 
tervals will essentially remain unchanged (will decrease 
by 1). Hence, in this case n oc T + . Similar arguments 
apply when T~ = T — T + <C T , leading to the proposed 
scaling. We neglect the second contribution, although it 
is not small. In Eq. 1(1911 we used n + oc (T + ) , while the 
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scaling part n + oc ((T — T + )/T) e was absorbed in the 
coefficients. We also should ideally recover the relation 



SH17T0 6 



which leads to 



sin tt9 

TT0 



■{*) 



£g(z)dz, 



(Dl) 



where z = T + /T and the factor of 2 arises because 
n + ~ n/2. In case of Eq. Ill9fl we have n + ~ 
(Tz) e sin7r6>/(7r6>) and relation {QUI is fulfilled approxi- 
mately. One can instead approximate n + ~ a(Tz) e or 
maybe n + ~ 6(1 — z) e (Tz) e where a or 6 will be deter- 
mined from Eq. IIDlll . Alternatively, a or 6 can be de- 
termined by equating the ensemble average of Eq. 1221) 
for small r with Eq. JT2J) for P+ = 1/2. 

There are two noticeable shortcomings of Eq. p2f . 
One of them regarding the discontinuity of its deriva- 
tive is mentioned in Section IVTl The other one becomes 
clear if one considers the complementary intensity signal 
J(t) = 1 - I(t). It follows from Eq. © that 



Ci(t',T') 



L [0,T] 



t',T'~ 



l + Cj(t',T') 



(D2) 



where Ci(t', T') = C TA (t', T') and Cj(t\ T') is the time- 
averaged correlation of signal J(t). Eq. Il22fl is written 
for Ci(t', T'), but analogous equation can be written for 
Cj(t' 7 T') as well, where T[o,t] is replaced by 1 — Zjo.T]- 
Then, unfortunately, the relation ilD2Jl will not hold in 
general. It will be satisfied trivially if t' = 0, or if t' 
is large enough so that one can use the second line of 
Eq. (EH for both Ci(t',T') and Cj(t',T') (more pre- 
cisely, if Ci(t', T") ~ Z[o T]1[t' T'l as m Eq. (12ofl and also 

Cj{t'X)^{i-i %T] ){i-i [tl>TI] )). 



Appendix E: BOUNDARIES OF C T A{t',T') 

Let us first consider the simpler case of r > 1/2: then 



T > t'. If T, 



[0,T'\ 



= V - T, 



[0,3" 



< t' — T, or equiva- 



lent^ I[o,T'] > 2(1 — r) then all the off intervals can 
be placed inside the interval [T, t'] and hence Cta^' , T') 
can attain its maximal value of 1, which we will write 
as Cta(P ,T') < 1, meaning that the limit is achievable. 
For Z[o,T'] < 2(1 — r) we put maximal duration of the 
off intervals inside the unused region [T, t'], making it 
identically zero, and the rest distribute identically on in- 
tervals [0,T] and \t',T'], so that all off intervals in [0,T] 



will be multiplied by all off intervals in [t',T r \, Then we 

have C TA {t',T') < (t - - (f - T)) /2) /T = 

X[ ^T']/(2(1— r )), where again, this upper bound is achiev- 
able. Considering the lower bound, for T]~ T ,j > T or 

equivalently 1[q,t'} < r , Cta(£ ' ,T') > and can reach 
zero, because we can make the whole interval [0,T] zero. 

For I [0 , T '] > r we have C T A{t', T') > (t — T^ Tl] ) /T = 
(2"[q jv| — r) /(l — r). Summarizing for r > 1/2 we have 
Eq.'£p. 

The case of r < 1/2, when t' < T, is more complicated. 
We note that if the interval lies inside of [t' ', T] it will be 
used twice, by both functions I(t) and I(t + i'). There- 
fore, to minimize CTA(t',T') for a given TjJ T ,j it seems 
desirable to put as much as possible of the off intervals 
into [i',? 1 ]. This is a good idea until we can make these 
intervals to be multiplied by the on intervals. If there 
is too much off time inside [t' ,T] then some zeros inside 
[t',T] will necessarily multiply other zeros inside 
the situation we want to avoid. This can happen only 
if r < 1/3. Therefore let us consider only the case of 
1/3 < r < 1/2. Then if T,~ < T - t' or equivalently 



[0,3"] 

l [0tT n > 2r yields C T A(t',T) > (T-2T, 



- [0,T'[ 



IT 



(2Z[ 0j t'] — T — l) /(l — r). For X[o.T'] — 2r we have 

C T A(t',T') > (t - 2(T - t') - (T [ - T , ] - (i - t'))) /T = 

(X[o,t'] + — l) /(l — r), assuming that if this bound 
is negative it is replaced by 0. For the upper 
bounds it follows that if TZ T n < 2[T — 2(T — t')] or 

I[o,T'] > 3 - 6r then C TA (t', T) < (t - T^ TI] /i) /T = 

(X [0 ,t'] + 1 - 2r) /(2(1 - r)) and if T+ Tr] < 3(T - f ) 

or I [0 , T /] < 3 - 6r then C TA {t' ,T) < (2T+ T ,,/3) /T = 

2T[q iT /]/(3(1 — r)). Summarizing for 1/3 < r < 1/2 yields 
Eq. '(Ell. 

Finally, for small r consider a simple counter-example 
showing that the bound C T A(t', T') < (l[o,T>] - r) /(l - 
r) can be overcome, in principle, for any 2jo,T']- Let 
T'/t' = 1/r be an integer. For any T^ T ,-, we then can 
distribute the on and o/f times by first filling the in- 
terval [0,f] with on time from to rTjJ T ,j and filling 

the remainder (from rT^ T ,^ to t' — rT 1 ) with off time. 
Rest of the intervals, [t',2t'], [2t',3t'], [T, T'] are filled 
in exactly the same way. Then clearly CTA(t',T r ) = 

(rT+ r/] (l - r)/r) /T = 1 [0 , T '] > (%t'] - r) /(l - r) 

for 2T[o,T'] < 1- The value 2"[o,t'] is n °t an upper bound 
either, in general, as can be seen, e.g., from Eq. Il-Tll) . 
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